function [y,ind,dist,yold,dens] = TSVQ(x,N,lowmem)
if nargin<3, lowmem=false; end
dmax = ceil(log2(N));
% Select the centroid of the entire dataset to be the root of the tree
y = mean(x);
yold = cell(1);
dens = size(x,1);d=0;
while size(y,1)<N
d = d+1;
% Perturb
y = splitTree(y,dmax,N);
if ~lowmem, yold{2*d-1} = y; end

iter=1;
while iter<10
iter = iter+1;
% Nearest centroids partition
mindist = ones(size(x,1),2)*Inf;
for i=1:size(y,1)
for j=1:size(x,1)
dist = sqrt(sum((x(j,:)-y(i,:)).^2));
if dist<mindist(j,1)
mindist(j,1) = dist;
mindist(j,2) = i;
end
end
end
% Centroid computation
dens = [];
for i=1:size(y,1)
dens(i) = sum(mindist(:,2)==i);
if dens(i)>0
y(i,:) = mean(x(mindist(:,2)==i,:));
end
end
%y = y(dens>0,:);
dens = dens(dens>0);
end
if ~lowmem, yold{2*d} = y; end
fprintf('.')
end
fprintf('\n')
ind = mindist(:,2);
dist = mindist(:,1);
function y = splitTree(y,dmax,N)
n = size(y,1);
vec = 1:n;
if N-n < (N/2)
vec = randperm(n);
vec = vec(1:(N-n));
end
for i = 1:length(vec)
y(n+i,:) = y(vec(i),:) + randn(1,size(y,2))/100;
end

